% 
%
%  function trial_time_list_ms = analyzeTiming( area_name, type_name,
%  phase_start_msec, phase_end_msec, phase_duration_msec , start_time_msec, end_time_msec)
% 
%  phase_start_msec, phase_end_msec, phase_duration_msec:   the start and
%		end times (inclusive) within the trial to analyze the data, and the phase_duration.  Start counting at 0 for the first msec.
%  start_time_msec, end_time_msec: start and end times of the experiments.
%
%  RETURNS:
%   A list of the time (in msec) at which the firing rate went above 2 Hz
%   in a 20 msec bin for any neuron in the group.
%   If threshold not exceeded, returns (phase_end_msec-phase_start_msec)+1;

%
%  EXAMPLE:
%
%
% analyzeTiming('S1','p23',200,799,1000,24000,40000)
%
%  Does not require a patternhist.dat file.

function trial_time_list_ms = analyzeSequence( area_name, type_name, phase_start_msec, phase_end_msec , trial_duration_msec, start_time_msec, end_time_msec)

read_groups; % Assumes groups.dat has the neural area data

index = [];
for i=1:length(areas)
    if strcmp( deblank(areas{i}), area_name) && strcmp( type_name, deblank(celltype{i})) 
        index = i;
        break;
    end
end
if isempty(index)
    disp([area_name ' ' type_name ' is not in groups.dat.']);
    return;
end

first = neuron_id(index,1); 
last  = neuron_id(index,2); 

num_neurons = max(neuron_id(:,2))+1;

first = first + 1;
last = last + 1;
rows = la(index);
cols = lb(index);
N=rows*cols;


save_seconds = 1;

% Move in short time windows and look for neural responses.
WINDOW = 20;  % msec.

filename = '';
% WINDOW should evenly divide 1000.

trial = 1;
for t1 =start_time_msec:trial_duration_msec:end_time_msec
    spikes = zeros(1, ceil((phase_end_msec-phase_start_msec)/WINDOW));
    i=1;

    for t = (t1+phase_start_msec):WINDOW:(t1+phase_end_msec)
        
        % find the file we need to open
        current_filename = ['spikes' num2str( ceil( (t+1)/(save_seconds*1000))*save_seconds*1000,15)  '.dat'];
        
        % Load new filename if necessary.
        if strcmp(current_filename, filename) == 0
            filename = current_filename;
            sdat = load(filename);
            max_time_msec = ceil( (t+1)/(save_seconds*1000))*save_seconds*1000;
            S=sparse( sdat(:,2)+1, sdat(:,1)+1, 1, max(num_neurons,max(sdat(:,1)+1)), max(max_time_msec,max(sdat(:,1))) );
        end
        
        % If this window is split across files, ignore it.
        end_filename = ['spikes' num2str( ceil( ( t+1 +WINDOW-1)/(save_seconds*1000))*save_seconds*1000,15)  '.dat'];
        if strcmp(end_filename, filename) == 0
            disp('warning: ignoring data because the bin window crossed file boundaries.');
            disp(t)
            disp(t+WINDOW-1)
            continue;
        end
        
        %disp('****************************');
        timebin = full(sum( S( first:last, t:(t+WINDOW-1)), 2));
        
        % We have the mfr in this time bin.  Now we have to see which training
        % pattern it matches.
        spikes(i) = sum(timebin);
        
        i=i+1;
    end
    
    indices = find(spikes);
    if ~isempty(indices)
        first_bin = indices(1);
        trial_time_list_ms(trial) = WINDOW*(first_bin-1);
    else
        trial_time_list_ms(trial) = (phase_end_msec-phase_start_msec)+1;
    end
    trial = trial + 1;
    
end

